library(stargazer)

# allows stargazer to use output from clm/clmm models
fake.clm <- function(x) {
  model.fake <- x
  model.fake$nobs <- x$dims$nobs
  model.fake$call[1] <- call("clm")
  model.fake
}

# estimate models
imm.model <- lmer(imm.scale ~ fb_lag*fb_d + gender + age + educ + faminc + pid7 + ideology + gender + own + hisp + asian + lived + pcthisp_01 + pctasian_01 + pctunemp_01 + medianhsvl_01 + (1|lookupzip), data = df.m, subset = imm == 1)

trump.model <- glmer(trump ~ fb_lag*fb_d + gender + age + educ + faminc + pid7 + ideology + gender + own + hisp + asian + lived + pcthisp_01 + pctasian_01 + pctunemp_01 + medianhsvl_01 + (1|lookupzip), data = df.m, subset = imm == 1, family = binomial(link="logit"), nAGQ = 0)

# produce tables
stargazer(imm.model, trump.model, out = "table4.html", star.cutoffs = c(.20, .10), digits = 2, notes.append = FALSE, single.row = TRUE)

# simulate QOIs for immigration opinion model
beta.hat <- fixef(imm.model)
vcov <- as.matrix(summary(imm.model)$vcov)
coefsim <- rmvnorm(1000, beta.hat, vcov)
k <- length(beta.hat)
X <- df.m[,c("gender","age","educ","faminc","pid7","ideology","own","hisp","asian","lived","pcthisp_01","pctasian_01","pctunemp_01","medianhsvl_01")]

fb.qs <- quantile(df.m$fb_d, c(0, 1), na.rm = TRUE)
lag.qs <- quantile(df.m$fb_lag, c(.05, .95), na.rm = TRUE)

rep.row<-function(x,n){
  matrix(rep(x,each=n),nrow=n)
}

means0 <- cbind(1, lag.qs[1], seq(fb.qs[1], 
                                  fb.qs[2], length.out = 1000), 
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[1] * seq(lag.qs[1], lag.qs[2], length.out = 1000))

results0 <- coefsim %*% t(means0)

means1 <- cbind(1, lag.qs[2], seq(fb.qs[1], 
                                  fb.qs[2], length.out = 1000),
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[2] * seq(lag.qs[1], lag.qs[2], length.out = 1000))

results1 <- coefsim %*% t(means1)

fd <- results1 - results0

z <- seq(fb.qs[1], fb.qs[2], length.out = 1000)

hist_values <- sample(z, 1000)

means <- apply(fd, 2, mean)
sdUpper <- apply(fd, 2, function(x) quantile(x, .95))
sdLower <- apply(fd, 2, function(x) quantile(x, .05))
z = seq(min(df.m$fb_d.n, na.rm = T), max(df.m$fb_d.n, na.rm = T), length.out = 1000)

hist_values <- c(0, diff(sapply(1:1000, function(x) as.numeric(sum(z[x] > df.m$fb_d.n, na.rm = T)))))

df <- data.frame(z = z, means = means * 100, 
                 sdUpper = sdUpper * 100, sdLower = sdLower * 100, hist_values) 

# generate plots
df %>% ggplot(aes(x = z, y = means)) + geom_line() + geom_ribbon(aes(ymin=sdLower,ymax=sdUpper),alpha=0.5) + geom_hline(yintercept = 0) + theme_bw() + labs(x = expression(Delta*" Immigrant (2000 - 2015)"), y = "First Difference") + geom_point(aes(z[1], means[1]), size = 10, shape = 21, fill = "white") + geom_point(aes(z[1000], means[1000]), size = 10, shape = 21, fill = "white") + geom_text(aes(z[1], means[1], label = round(means[1], 0), family = "Palatino Linotype")) + geom_text(aes(z[1000], means[1000], label = round(means[1000], 0), family = "Palatino Linotype")) + theme(plot.title = element_text(size=13, family="Palatino Linotype", face = "bold"), axis.title = element_text(size = 12, family = "Palatino Linotype")) + geom_density(aes(x = z, y = hist_values/100), stat = "identity", fill = "black")


ggsave("pb_fight2.pdf", device=cairo_pdf, dpi = "retina")

# simulate QOIs for Trump model
beta.hat <- fixef(trump.model)
vcov <- as.matrix(summary(trump.model)$vcov)
coefsim <- rmvnorm(1000, beta.hat, vcov)
k <- length(beta.hat)
X <- df.m[,c("gender","age","educ","faminc","pid7","ideology","own","hisp","asian","lived","pcthisp_01","pctasian_01","pctunemp_01","medianhsvl_01")]

fb.qs <- quantile(df.m$fb_d, c(0, 1), na.rm = TRUE)
lag.qs <- quantile(df.m$fb_lag, c(.05, .95), na.rm = TRUE)

rep.row<-function(x,n){
  matrix(rep(x,each=n),nrow=n)
}

means0 <- cbind(1, lag.qs[1], seq(fb.qs[1], 
                                  fb.qs[2], length.out = 1000), 
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[1] * seq(lag.qs[1], lag.qs[2], length.out = 1000))

results0 <- coefsim %*% t(means0)
results0 <- apply(results0, 1:2, function(x) invlogit(x))

means1 <- cbind(1, lag.qs[2], seq(fb.qs[1], 
                                  fb.qs[2], length.out = 1000),
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[2] * seq(lag.qs[1], lag.qs[2], length.out = 1000))

results1 <- coefsim %*% t(means1)
results1 <- apply(results1, 1:2, function(x) invlogit(x))

fd <- results1 - results0

z = seq(min(df.m$fb_d.n, na.rm = T), max(df.m$fb_d.n, na.rm = T), length.out = 1000)

means <- apply(fd, 2, mean)
sdUpper <- apply(fd, 2, function(x) quantile(x, .95))
sdLower <- apply(fd, 2, function(x) quantile(x, .05))

df <- data.frame(z = z, means = means * 100, 
                 sdUpper = sdUpper * 100, sdLower = sdLower * 100) 
# generate plots
df %>% ggplot(aes(x = z, y = means)) + geom_line() + geom_ribbon(aes(ymin=sdLower,ymax=sdUpper),alpha=0.5) + geom_hline(yintercept = 0) + theme_bw() + labs(x = expression(Delta*" Immigrant (2000 - 2015)"), y = "First Difference") + geom_point(aes(z[1], means[1]), size = 10, shape = 21, fill = "white") + geom_point(aes(z[1000], means[1000]), size = 10, shape = 21, fill = "white") + geom_text(aes(z[1], means[1], label = round(means[1], 0), family = "Palatino Linotype")) + geom_text(aes(z[1000], means[1000], label = round(means[1000], 0), family = "Palatino Linotype")) + theme(plot.title = element_text(size=13, family="Palatino Linotype", face = "bold"), axis.title = element_text(size = 12, family = "Palatino Linotype")) + geom_density(aes(x = z, y = hist_values/50), stat = "identity", fill = "black")

ggsave("pb_fight3.pdf", device=cairo_pdf, dpi = "retina")

